arXiv:l 608.0068 lvl  [quant-ph]  2  Aug  2016 


Observation  of  Prethermalization  in  Long-Range  Interacting  Spin  Chains 


B.  Neyenhuis,1^  J.  Smith,1  A.  C.  Lee,1  J.  Zhang,1  P.  Richerme,2 
P.  W.  Hess,1  Z.-X.  Gong,1  A.  V.  Gorshkov,1  and  C.  Monroe1 

1  Joint  Quantum  Institute  and  Joint  Center  for  Quantum  Information  and  Computer  Science, 
University  of  Maryland  Department  of  Physics  and  National  Institute 
of  Standards  and  Technology,  College  Park,  Maryland  20742 
2 Department  of  Physics,  Indiana  University,  Bloominqton,  IN,  47405 
(Dated:  August  3,  2016) 

Statistical  mechanics  can  predict  thermal  equilibrium  states  for  most  classical  systems,  but  for 
an  isolated  quantum  system  there  is  no  general  understanding  on  how  equilibrium  states  dynam¬ 
ically  emerge  from  the  microscopic  Hamiltonian  m-  For  instance,  quantum  systems  that  are 
near-integrable  usually  fail  to  thermalize  in  an  experimentally  realistic  time  scale  and,  instead,  relax 
to  quasi-stationary  prethermal  states  that  can  be  described  by  statistical  mechanics  when  approx¬ 
imately  conserved  quantities  are  appropriately  included  in  a  generalized  Gibbs  ensemble  (GGE) 
EM!!-  Here  we  experimentally  study  the  relaxation  dynamics  of  a  chain  of  up  to  22  spins  evolving 
under  a  long-range  transverse  field  Ising  Hamiltonian  following  a  sudden  quench.  For  sufficiently 
long-ranged  interactions  the  system  relaxes  to  a  new  type  of  prethermal  state  that  retains  a  strong 
memory  of  the  initial  conditions.  In  this  case,  the  prethermal  state  cannot  be  described  by  a  GGE, 
but  rather  arises  from  an  emergent  double-well  potential  felt  by  the  spin  excitations.  This  result 
shows  that  prethermalization  occurs  in  a  significantly  broader  context  than  previously  thought, 
and  reveals  new  challenges  for  a  generic  understanding  of  the  tliermalization  of  quantum  systems, 
particularly  in  the  presence  of  long-range  interactions  El- 


Statistical  mechanics  can  predict  thermal  equilibrium 
states  for  most  classical  systems,  but  for  an  isolated  quan¬ 
tum  system  there  is  no  general  understanding  on  how 
equilibrium  states  dynamically  emerge  from  the  micro¬ 
scopic  Hamiltonian.  For  instance,  quantum  systems  that 
are  near-integrable  usually  fail  to  thermalize  in  an  ex¬ 
perimentally  realistic  time  scale  and  instead,  relax  to 
quasi-stationary  prethermal  states  that  can  be  described 
by  statistical  mechanics  when  approximately  conserved 
quantities  are  appropriately  included  in  a  generalized 
Gibbs  ensemble  (GGE).  Here  we  experimentally  study 
the  relaxation  dynamics  of  a  chain  of  up  to  22  spins  evolv¬ 
ing  under  a  long-range  transverse  field  Ising  Hamiltonian 
following  a  sudden  quench. 

In  the  classical  world  thermalization  is  expected  in 
all  but  special  cases  where  conserved  quantities  or  hid¬ 
den  symmetries  prevent  the  ergodic  exploration  of  phase 
space.  Because  the  classical  world  is  ultimately  com¬ 
prised  of  quantum  systems  we  therefore  expect  that  a 
closed  quantum  system  will  also  reach  thermal  equilib¬ 
rium.  Although  quantum  dynamics  are  unitary,  mea¬ 
surements  made  within  a  subsystem  trace  over  the  rest 
of  the  system  and  appear  thermal  because  the  rest  of  the 
system  acts  as  a  thermal  bath  mm- 

However  this  is  not  always  the  case.  For  integrable 
models  an  extensive  number  of  conserved  quantities  pre¬ 
vent  the  efficient  exploration  of  phase  space  SSI  and 
the  system  relaxes  to  a  steady-state  predicted  by  a  GGE 
mm  specified  by  the  initial  values  of  the  integrals  of 
motion.  For  near-integrable  systems,  such  as  weakly  in¬ 
teracting  ultracold  gases,  thermalization  can  still  occur, 
but  only  over  extremely  long  time  scales  beyond  current 
experimental  reach  mm-  However,  it  is  possible  to 


observe  quasi-stationary  states,  often  called  prethermal, 
that  emerge  within  an  experimentally  accessible  time 
scale. 

Previous  observations  of  prethermal  states  have  fo¬ 
cused  on  those  described  by  a  GGE  associated  with  the 
integrable  part  of  the  model  (10;,  111.  Here,  we  observe  a 
new  form  of  prethermalization  m,  where  a  system  of  in¬ 
teracting  spins  rapidly  evolves  to  a  quasi-stationary  state 
that  cannot  be  predicted  by  a  GGE.  This  type  of  prether¬ 
mal  state  arises,  even  in  the  thermodynamic  limit,  when 
a  system  has  long-range  interactions  and  open  bound¬ 
aries  such  that  the  translational  invariance  is  broken.  As 
a  result,  spin  excitations  feel  an  emergent  double- well  po¬ 
tential  whose  depth  grows  with  interaction  range  (Fig|T|) . 
Memory  of  the  initial  state  is  preserved  by  this  emergent 
potential,  but  is  eventually  lost  due  to  weak  tunneling 
between  the  two  wells  and  the  interactions  between  spin 
excitations. 

Effective  spin- 1/2  particles  are  encoded  in  the 
2Si/2|E  =  0,m_p  =  0)  and  |F  =  1,to_f  =  0)  hyperfine 
‘clock’  states  of  a  mYb+  ion,  denoted  |!)2  and  |f)z  [20]- 
We  confine  a  chain  of  ions  in  a  linear  rf  Paul  trap  and 
apply  optical  dipole  forces  to  generate  the  effective  spin- 
spin  coupling  mm  of  an  Ising  Hamiltonian  (methods) : 

h  =  c1) 

i<j  i 

where  ex/  (7  =  x ,  z)  are  the  Pauli  matrices  acting  on 
the  ith  spin,  is  the  coupling  between  spins  i  and 
j,  B/(2tt)  =  10  kHz  is  a  uniform  effective  transverse 
field,  and  we  use  units  in  which  Plancks  constant  equals 
1  (methods).  The  spin-spin  interaction  is  long-range  and 
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can  be  described  by  a  power  law  decay  where  Jy  = 
\i-j\c  i  dmax/( 27t)  is  the  maximum  coupling  strength 
which  ranges  from  0.45  kHz  to  0.98  kHz.  We  tune  the 
power  law  exponent  a  between  0.55  and  1.33  by  chang¬ 
ing  the  axial  confinement  of  the  ions.  With  long-range 
interactions,  H  is  in  general  non-integrable  (in  contrast 
to  the  nearest-neighbor  case  where  the  ID  model  is  inte- 
grable  |23j)  and  thermalization  is  anticipated  in  the  long 
time  limit  according  to  the  eigenstate  thermalization  hy¬ 
pothesis  fMWM .  However,  one  can  map  (jTJ)  with  a  small 
number  of  spin  excitations  into  a  near-integrable  model 
of  bosons,  with  an  integrable  part  made  of  free  bosons 
and  an  integrability-breaking  part  consist  of  weak  inter¬ 
actions  among  bosons  that  leads  to  the  eventual  thermal¬ 
ization. 

We  initialize  the  chain  of  seven  spins  by  optically 
pumping  all  spins  to  the  state,  and  then  use  a  tightly 
focused  individual-ion  addressing  laser  to  excite  a  single 
spin  on  one  end  of  the  chain  to  the  |f)-  state  as  seen  in 
Fig.  [T}i  (methods)  [27].  The  spins  then  evolve  under  0 
and  we  measure  the  time  evolution  of  the  spin  projec¬ 
tion  in  the  z-basis.  For  the  shortest  range  interactions 
we  realize  ( a  =  1.33)  the  system  rapidly  evolves  to  a 
prethermal  state  predicted  by  the  GGE  associated  with 
the  integrals  of  motion  corresponding  to  the  momentum 
space  occupation  number  of  the  non-interacting  bosons, 
which  does  not  preserve  memory  of  the  initial  spin  exci¬ 
tation  location  (methods)  (Fig  lb). 

However,  in  the  long-range  interacting  case  ( a  =  0.55) 
we  see  the  position  of  the  spin  excitation  reaches  an  equi¬ 
librium  value  that  retains  a  memory  of  the  initial  state 
(Fig.  §1)  out  to  the  longest  experimentally  achievable 
time  of  25/Jmax.  This  prethermal  state  is  in  obvious 
disagreement  with  both  a  thermal  state  and  the  GGE 
prediction,  which  both  maintain  the  right-left  symmetry 
of  the  system. 

The  dynamics  of  the  spin-wave  boson  model  for  short- 
range  interactions  are  similar  to  those  of  a  free-particle 
in  a  square-well  potential.  However,  long-range  interac¬ 
tions  distort  the  square-well  to  a  double-well  potential 
Fig.  [l}r.  Here  we  emphasize  that  the  double-well  poten¬ 
tial  emerges  non-trivially  because  our  model  ([I])  is  tran¬ 
sitionally  invariant  without  boundaries.  The  spin-wave 
boson  model  has  an  extensive  number  of  near-degenerate 
eigenstates  that  are  symmetric  and  antisymmetric  super¬ 
positions  of  spin  excitations  in  the  left  and  right  po¬ 
tential  wells.  For  seven  spins,  we  calculate  the  energy 
difference,  A£mn,  between  all  pairs  of  eigenstates  and 
plot  them  with  respect  to  a  in  F ig jTJ-  where  the  ampli¬ 
tudes,  \{m\ipo)  (ipo\n)\2 ,  are  products  of  the  overlaps  ofthe 
eigenstates  |m)  and  |?r)  overlap  with  the  initial  statej^o)- 
With  a  =  0.55  the  two  lowest  energy  states  are  almost 
degenerate,  with  an  energy  difference  of  approximately 
one-thousand  times  smaller  than  Jmax,  due  to  the  tun¬ 
neling  rate  between  the  ground  states  of  each  well,  which 
is  exponentially  small  in  the  barrier  height.  As  a  result 


FIG.  1.  Emergent  double  well  potential,  (a)  The  spin 
chain  starts  with  a  single  spin  excitation  on  the  left  end  in 
an  effective  double-well  potential,  Ueff,  whose  barrier  height 
is  determined  by  the  range  a  of  the  interactions,  (b)  For 
short  range  interactions  (a  =  1.33),  we  map  the  system  to 
a  particle  in  a  ID  square  well  where  the  excitation  becomes 
symmetrically  distributed  across  the  chain  as  predicted  by 
the  GGE,  {cf) gge •  (d)  However,  for  long-range  interactions 
(a  =  0.55)  there  is  an  emergent  double-well  potential  which 
prevents  the  efficient  transfer  of  the  spin,  and  the  excitation 
location  retains  memory  of  the  initial  state,  in  contrast  to 
( oI)gge-  (c)  The  double  well  gives  rise  to  near-degenerate 
eigenstates  as  a  is  decreased  as  seen  in  the  calculated  energy 
difference  between  all  pairs  of  eigenstates  versus  a,  with  am¬ 
plitude  weighted  by  the  product  of  the  eigenstates’  overlap 
with  the  initial  state. 


the  spin  excitation  will  remain  in  its  initial  well  until  it 
tunnels  across  the  potential  barrier  at  much  longer  times. 

To  better  characterize  the  “location”  of  the  spin  exci¬ 
tation  we  construct  the  observable 


__  ^  2i  —  N  -  1  a-  +  1 
^  N-  1  2 


(2) 


where  N  is  the  number  of  ions.  The  expectation  value 
of  C  varies  between  -1  and  1  for  a  spin  excitation  on 
the  left  and  right  ends,  respectively.  Due  to  the  spatial 
inversion  symmetry  of  the  spin-wave  boson  model  and 
0.  both  the  GGE  predicted  and  thermal  values  of  ( C ) 
are  zero.  In  Fig.  [2]we  plot  (C)  along  with  its  cumulative 
time  average  ( C )  for  a  =  1.33  (short-range)  and  a  =  0.55 
(long-range).  To  further  accentuate  the  asymmetry  of  the 
prethermalization  we  prepare  initial  states  with  a  single¬ 
spin  excitation  on  the  left  or  right  ends  of  the  chain. 

With  short-range  interactions  the  system  quickly  re¬ 
laxes  to  the  prethermal  value  predicted  by  the  GGE 
where  ( C )  is  near  zero  irrespective  of  the  initial  state. 
But  with  long-range  interactions,  the  prethermal  state 
that  emerges  retains  a  clear  memory  of  the  initial  con¬ 
ditions  and  is  different  than  both  the  thermal  and  GGE 
predictions. 

It  is  useful  to  talk  about  thermalization  in  terms  of  lo¬ 
cal  quantities  since  subsystems  must  use  the  rest  of  the 
chain  as  a  heat  bath.  In  Fig.  [2]  we  plot  the  cumulative 
time  average  of  the  deviation  of  the  single  spin  magne- 
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FIG.  2.  Measured  location  of  spin  excitation.  The  average  position  of  the  spin  excitation,  (C),  is  plotted  for  an  initial 
excitation  on  the  left  (dark  blue  unfilled  squares)  and  right  (dark  red  unfilled  circles)  along  with  their  cumulative  time  average 
(blue  and  red  filled  circles  and  squares),  for  short-range  (a  =  1.33)  and  long-range  ( a  =  0.55)  interactions  with  initial  states 
with  one  spin  excitation  (top  panels)  and  two  spin  excitations  (bottom  panels).  The  cumulative  time  average  of  the  deviation 
of  the  postselected  individual  spin  projections  from  the  generalized  Gibbs  ensemble,  (<rf )  —  (af)GGE,  is  also  plotted.  For  short- 
range  interactions  the  spins  quickly  thermalize  to  a  value  predicted  by  the  GGE,  but  for  the  long-range  interacting  system  a 
new  type  of  prethermal  state  emerges.  Error  bars,  1  s.d. 


tizations,  {erf),  from  the  equilibrium  value  predicted  by 
the  GGE.  Here  we  postselect  the  data  for  the  correct 
number  of  spin  excitations  to  eliminate  the  effects  of  im¬ 
perfect  state  preparation  and  small  deviations  from  our 
model  Hamiltonian  due  to  unwanted  excitations  of  the 
phonon  modes  (methods).  For  the  short-range  interac¬ 
tions  we  see  the  cumulative  time  average  for  each  spin 
quickly  converge  to  the  GGE  predicted  value.  However, 
with  the  long-range  interactions  we  see  the  individual 
spins  reach  a  steady  state  that  does  not  match  the  GGE 
as  the  emergent  double- well  prevents  the  efficient  transfer 
of  the  excitation  across  the  chain. 

We  show  the  prethermal  state’s  robustness  to  weak  in¬ 
teractions,  similar  to  many  body  localization  0E],  by 
preparing  initial  states  with  two  spin  excitations.  In  this 
case,  the  multiple  spin  flips  increase  the  size  of  the  in- 
tegrability  breaking  part  of  the  Hamiltonian  which  rep¬ 
resents  weak  interactions  between  the  spin-wave  bosons. 
However,  the  prethermal  state  still  persists.  For  bet¬ 
ter  contrast  between  the  prethermal  and  GGE  predicted 
values  of  (G),  we  flip  the  second  and  fourth  spin  such 
that  |i/>o)  =  |4-t4't'l'44') •  We  emphasize  that  the  observed 
prethermalization  and  departure  from  the  GGE  predic¬ 
tion  is  not  sensitive  to  the  specific  choice  of  initial  state  in 
the  thermodynamic  limit.  But  for  a  small  sized  system, 
these  initial  states  offer  us  the  maximum  signal.  As  be¬ 
fore  we  also  prepare  the  mirror  image  of  the  initial  state 
by  exciting  the  fourth  and  sixth  ions  (IV’o)  =  |4~l4-t4-'N'))- 
We  observe  relaxation  to  the  value  predicted  by  the  GGE 
for  short-range  interactions,  but  with  long-range  interac¬ 
tions  we  see  a  prethermal  state  that  strongly  deviates 
from  the  GGE  (bottom  panel  of  Fig.  [2]). 


FIG.  3.  General  feature  of  the  model.  Comparison  be¬ 
tween  a  numerical  simulation  of  the  transverse  field  Ising 
model  and  experimental  data  for  single  (a)  and  double  (b) 
initial  spin  excitations.  The  inset  shows  excellent  agreement 
between  experiment  (black  dots)  and  numerics  (red  line)  ac¬ 
counting  for  experimental  noise  for  (C).  For  both  the  single 
and  double  spin  flip  cases,  the  prethermal  state  persists  for 
much  longer  than  the  experimental  timescale  and  eventually 
relaxes  to  the  thermal  state.  Error  bars,  1  s.d. 


In  Fig.  [3]  we  plot  the  experimental  evolution  of  the 
prethermal  state  for  both  the  double  and  single  spin  flip 
initial  states  along  with  a  long-time  numerical  simula¬ 
tion  under  ([!])  accounting  for  known  experimental  noise 
(methods).  The  plots  demonstrate  excellent  agreement 
between  numerical  simulations  and  experimental  data 
and  confirm  that  the  prethermal  states  persist  well  be¬ 
yond  the  current  experimental  time  limit.  Due  to  the 
non-conservation  of  the  number  of  spin-wave  bosons  for 


4 


1.0 

22  spins  all  up 

0.5 

Single  spin  excitation 

Q°.° 

-0.5 

i/>- -  f 

End  of  spin  evolution  / 

-1.0 

0  6  12  18  24  30  36 

Jmax  t 

FIG.  4.  Scaling  to  larger  system  size.  Time  evolution 
(light  blue)  and  cumulative  time  average  (orange)  of  (C)  with 
false  color  pictures  of  the  22  ion  chain,  where  the  brightness  of 
each  ion  is  determined  by  the  value  of  (of).  The  ions  fluoresce 
during  detection  when  in  the  |f)2  state  (top  picture).  We 
initialize  the  spins  with  a  single  spin  excitation  on  the  left 
end  (middle  picture).  After  evolving  for  36  Jmax  the  spin 
excitation  is  delocalized,  but  its  average  position  remains  on 
the  left  half  of  the  chain  (bottom  picture).  Error  bars,  1  s.d. 

any  finite  B  and  the  interactions  between  them,  the  sys¬ 
tem  will  eventually  relax  to  the  thermal  equilibrium  in 
the  thermodynamic  limit,  however  relaxation  to  the  GGE 
may  or  may  not  be  seen  depending  on  the  range  of  inter¬ 
actions  (methods). 

To  demonstrate  that  the  prethermal  state  we  observe  is 
not  sensitive  to  system  size,  we  repeat  the  experiment  in 
a  chain  of  22  ions,  the  largest  ion  chain  used  for  quantum 
simulation  in  the  literature  to  date.  The  time  evolution 
and  cumulative  time  average  of  ( C )  are  depicted  in  Fig. 
[4j  Experimentally,  as  well  as  in  the  analytic  result  in 
the  thermodynamic  limit  (methods),  we  see  the  system 
relaxes  to  a  similar  quasi-equilibrium  state  as  before  that 
clearly  has  memory  of  the  initial  state. 

We  point  out  the  observed  prethermalization  and  de¬ 
viation  from  the  GGE  should  disappear  if  the  system 
is  subject  to  periodic  boundary  conditions,  regardless  of 
system  size.  Here  we  emphasize  that  the  long-range  inter¬ 
actions  make  the  boundary  conditions  relevant  for  bulk 
properties,  and  thus,  changing  the  boundary  conditions 
can  impact  an  extensive  number  of  eigenstates.  The  ef¬ 
fects  of  long-range  interactions  on  many-body  dynamics 
are  far  richer  than  the  effect  discussed  in  the  current  ex¬ 
periment.  For  sufficiently  long-range  interactions,  the 
notion  of  locality  breaks  down  and  quasi-particles  in  the 
system  can  travel  at  divergent  velocities  for  thermody¬ 
namic  systems,  potentially  leading  to  dramatically  differ¬ 
ent  thermalization/prethermalization  time  scales  in  cer¬ 
tain  systems  p8H30l.  We  believe  that  the  current  exper¬ 
iment,  as  well  as  the  platform  it  is  built  upon,  will  pave 
the  way  to  a  more  complete  understanding  of  the  fun¬ 
damental  role  long-range  interactions  play  in  the  quench 


dynamics  and  emergent  statistical  physics  of  quantum 
many-body  systems. 
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METHODS 


Effective  Hamiltonian  Generation 

We  generate  spin-spin  interactions  by  applying  spin- 
dependent  optical  dipole  forces  to  ions  confined  in  a  3- 
layer  linear  Paul  trap  with  a  4.8  MHz  radial  frequency. 
Two  off-resonant  laser  beams  with  a  wavevector  differ¬ 
ence  5  k  along  a  principal  axis  of  transverse  motion  glob¬ 
ally  address  the  ions  and  drive  stimulated  Raman  transi¬ 
tions.  The  two  beams  contain  a  pair  of  beatnote  frequen¬ 
cies  symmetrically  detuned  from  the  resonant  transition 
at  r'o  =  12.642819  GHz  by  a  frequency  /i,  comparable  to 
the  transverse  motional  mode  frequencies.  In  the  Lamb- 
Dicke  regime,  this  results  in  the  Ising-type  Hamiltonian 
in  Eq.  (1)  [21]  |22|  E]  with 

j  _  h(5k)2n2  ^  vi>mvjtm 

lJ~  2  M  2-<’  W 

m= 1  m 

where  Q  is  the  global  Rabi  frequency,  H  is  the  reduced 
Planck’s  constant,  V)j?Ti  is  the  normal  mode  matrix  TT2] . 


and  (jjm  are  the  transverse  mode  frequencies.  The  cou¬ 
pling  profile  may  be  approximated  as  a  power-law  decay 
Jij  ss  Jmax/\i  —  j\a,  where  in  principle  a  can  be  tuned 
between  0  and  3  by  varying  the  laser  detuning  /i  or  the 
bandwidth  of  u>m  EUSEJ.  For  the  7  ion  data  in  this  work, 
a  was  tuned  to  0.55  (long-range  interactions)  and  1.33 
(short-range  interactions)  by  changing  the  bandwidth  of 
uim.  By  asymmetrically  adjusting  the  laser  beatnote  de¬ 
tuning  /r  about  the  carrier  by  a  value  of  2 B  we  apply  a 
uniform  effective  transverse  magnetic  field  of  Baf  [33]. 

Single  spin  flip  initialization  and  site-resolved 
detection 

We  initialize  individual  spin  excitations  using  a  tightly 
focused  laser  beam  to  imprint  a  fourth  order  AC  Stark 
shift  [27]  in  conjunction  with  a  Ramsey[34]  or  Rabi  se¬ 
quence  as  seen  in  Extended  Data  Fig.  ([!]).  When  the  ion 
spacing  is  larger  than  the  beam  waist  of  the  individual- 
ion  addressing  laser  as  is  the  case  for  the  7  ion  short  range 
data  we  employ  a  Ramsey  method.  This  consists  of  first 
optically  pumping  the  spins  to  |].-).  Then  we  globally 
perform  a  7r/2  rotation  so  that  all  of  the  spins  are  in 
|j,x).  Using  the  individual-ion  addressing  beam,  a  Stark 
shift  is  applied  to  the  spins  to  be  flipped,  and  then  we 
allow  the  chain  to  evolve  until  these  spins  are  tt  out  of 
phase  compared  to  the  spins  without  an  applied  Stark 
shift.  Afterwards,  a  global  7r/2  rotation  brings  the  spins 
back  into  the  z-basis.  With  this  method,  individual  spin 
flips  can  be  prepared  with  a  fidelity  of  ~  0.97,  while  N 
spin  flips  can  be  achieved  with  a  fidelity  of  ~  (0.97)^. 

We  employ  the  Rabi  method  for  the  long  range  in¬ 
teracting  data  because  site-resolved  Stark  shifts  can  no 
longer  be  applied  since  the  ion  separation  is  smaller  than 
the  beam  waist.  Here,  we  apply  a  large  Stark  shift  to 
all  of  the  spins  except  the  ones  to  be  flipped  and  apply 
a  global  7r  pulse  at  the  hyperfine  splitting  between  the 
two  effective  spin  levels.  Thus,  only  the  ions  without 
an  applied  Stark  shift  are  flipped.  This  approach  has  a 
single  and  N  spin  flip  fidelity  of  ~  0.85  and  ~  (0.85)^ 
respectively. 

After  quenching  to  and  allowing  time  evolution  under 
our  spin  Hamiltonian,  we  measure  the  spin  projections 
of  each  ion  along  the  z  direction  of  the  Bloch  sphere. 
We  expose  the  ions  to  a  laser  beam  that  addresses  the 
cycling  transition  2Si/2 1 W  =  1)  to  2Pi/2|E  =  0)  for  3  ms. 
Ions  fluoresce  only  if  they  are  in  the  state  |f)z-  This 
fluorescence  is  collected  through  an  NA=0.23  objective 
and  imaged  using  an  intensified  CCD  camera  with  single¬ 
site  resolution. 

To  discriminate  between  ‘bright’  and  ‘dark’  states  (|f)z 
and  |j,)_  respectively),  we  begin  by  calibrating  the  camera 
with  1000  cycles  each  of  all-bright  and  all-dark  states. 
For  the  bright  states,  the  projection  of  the  2D  CCD  image 
onto  a  one-dimensional  row  gives  a  profile  comprised  of 
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Extended  Data  Figure  1.  Top  box:  We  use  a  Ramsey  method 
of  preparing  individual  spin  flips  when  the  ion  spacing  is  larger 
than  the  individual  addressing  beam  waist.  The  spins  are 
first  optically  pumped  to  |4-)„,  and  then  a  global  7t/2  rotation 
brings  them  into  the  aty-plane  of  the  Bloch  sphere.  The  sys¬ 
tem  then  evolves  with  a  Stark  shift  applied  to  the  ions  to  be 
flipped  until  the  Stark  shifted  spins  are  7r  out  of  phase  com¬ 
pared  to  the  other  spins.  Then  a  global  7t/2  rotation  brings 
the  spins  back  into  the  2-basis.  Bottom  box:  We  use  a  Rabi 
method  of  preparing  spin  excitations  when  the  ion  spacing  is 
smaller  than  the  individual  addressing  beam  waist.  In  this 
case  the  ions  are  optically  pumped  to  |4.)z  and  then  a  large 
Stark  shift  is  applied  to  all  of  the  ions  except  ones  to  be  flipped 
while  a  global  7r  rotation  is  simultaneously  performed  at  the 
hyperfine  splitting  between  the  effective  spin  states. 


Gaussian  distributions  at  each  ion  location.  We  perform 
fits  to  locate  the  center  and  fluorescence  width  of  each 
ion. 

We  achieve  single-shot  discrimination  of  individual  ion 
states  in  the  experimental  data  by  fitting  the  captured 
one-dimensional  profile  to  a  series  of  Gaussian  distribu¬ 
tions  with  calibrated  widths  and  positions  but  freely- 
varying  amplitudes.  These  extracted  values  for  each  ion 
are  then  compared  with  a  threshold  found  via  Monte- 
Carlo  simulation  to  determine  whether  the  measured 
state  was  ‘bright’  or  ‘dark’.  Our  discrimination  protocol 
also  gives  an  estimate  of  the  detection  error  (e.g.  mis¬ 
diagnosing  a  ‘bright’  ion  as  ‘dark’),  which  is  typically  of 
order  ~  5%.  Corrected  state  probabilities  (along  with 
their  respective  errors)  are  found  following  the  method 
outlined  in  |.'15j ,  which  also  takes  into  account  errors  due 
to  quantum  projection  noise. 


Experimental  noise  sources  and  their  influence  on 
the  thermalization  dynamics. 

As  discussed  in  the  text,  there  are  fluctuations  on  the 
interaction  strength  Jjj,  which  originate  from  noise  on 
the  laser  intensity  and  wm  [3Hj  -  This  noise  is  slow  com¬ 


pared  to  a  single  experiment,  but  fast  compared  to  the 
thousands  of  experiments  it  takes  to  complete  a  data  set. 
Averaging  over  this  classical  noise  leads  to  dynamics  that 
resemble  a  running  time  average,  because  the  fast  tem¬ 
poral  oscillations  are  effectively  canceled  out  by  the  fluc¬ 
tuating  Jij.  To  account  for  this,  the  numeric  simulations 
average  over  a  small  range  of  coupling  strengths  (stan¬ 
dard  deviation  of  0.12  Jmax)- 

Another  source  of  noise  is  the  fourth  order  AC  Stark 
shift  from  the  Mplmer-Sprensen  interaction  laser  side¬ 
bands.  This  noise  also  has  a  negligible  effect,  since  we  are 
in  the  regime  of  large  transverse  fields  and  this  Stark  shift 
term  only  adds  a  small  global  az  fluctuation  of  about  30 
Hz,  on  top  of  the  10  kHz  transverse  field  applied.  We  ex¬ 
perimentally  verify  this  by  applying  a  global  offset  field 
of  ±  400  Hz  and  observe  no  difference  in  the  observed 
prethermalization.  The  relaxation  dynamics  are  robust 
against  these  experimental  noise  sources,  however  it  is 
sensitive  to  asymmetries  of  the  spin-spin  coupling  ma¬ 
trix. 


Measuring  the  spin-spin  coupling  matrix 

We  directly  measure  the  spin-spin  coupling  matrices 
for  seven  ions  for  both  long  and  short  range  interactions 
and  ensure  it  is  symmetric  as  seen  in  Extended  Data 
Fig.  §  In  order  to  measure  the  strength  of  the  inter¬ 
action  between  two  spins,  we  shelve  all  but  the  two  ions 
of  interest  out  of  the  interaction  space  and  directly  ob¬ 
serve  their  time  evolution.  This  is  done  by  first  perform¬ 
ing  individual  rotations  on  these  two  ions  to  the  |f)2  as 
outlined  above.  Then  we  perform  a  global  7r  rotation  be¬ 
tween  |4-)z,  2S1/2|F  =  0,  nip  =  0),  and  one  of  the  Zeeman 
states,  2S1/2|E  =  1,?jif  =  —1),  which  takes  the  other  5 
spins  out  of  the  interaction  space.  We  then  apply  the 
Hamiltonian  which  now  only  acts  on  the  two  remaining 
spins. 


Justification  for  postselection 

As  noted  above  the  initial  N  spin  flip  fidelity  is  ap¬ 
proximately  (0.97)^  and  (0.85)^  for  short  and  long  range 
interactions  respectively.  Numerically,  we  find  that  the 
number  of  spin  excitations  is  essentially  constant  on  the 
experimental  timescale  which  is  expected  because  the 
transverse  field  Ising  model  can  be  mapped  to  an  XY 
model  for  sufficiently  large  transverse  field  and  the  XY 
model  conserves  the  number  of  spin  excitations  251. 
However,  experimentally  we  observe  significant  leakage 
out  of  the  N  excitation  subspace  with  less  than  50%  re¬ 
maining  at  the  end  of  the  evolution.  Thus,  we  postse¬ 
lect  the  data  for  the  correct  number  of  spin  excitations 
to  eliminate  the  effects  of  imperfect  state  preparation, 
detection  error,  and  small  deviations  from  our  model 
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Extended  Data  Figure  2.  We  directly  measure  the  spin-spin 
coupling  matrix  with  7  ions  for  both  short  (left  matrix)  and 
long  (right  matrix)  range  interactions  and  see  it  is  symmetric. 
To  measure  the  coupling  between  two  spins  we  transfer  all 
except  the  spins  of  interest  out  of  the  interaction  space.  This 
is  done  by  shelving  the  other  spins  in  one  of  the  Zeeman  levels. 


Hamiltonian  due  to  unwanted  excitations  of  the  phonon 
modes. 


The  spin-boson  mapping  and  the  generalized  Gibbs 
ensemble 

To  explain  our  observed  prethermalization,  it  is  con¬ 
venient  to  map  the  spins  into  bosons  by  using  the 
Holstein-Primakoff  transformation:  of  =  2ajai  —  1,  of  = 

a\  ijl  —  a\  di  .  We  will  assume  that  the  average  spin  ex¬ 
citation  density  n  =  (>A(a|aj)/iV  is  much  smaller  than  1. 
This  assumption  is  justified  in  our  experiments  because 
our  initial  states  have  small  spin  excitation  densities  and 
we  set  max(Jij)  C  5  so  the  amount  of  n  that  will  be 
dynamically  created  is  small  [~  (ma x{Jjj)  /  B)2].  There¬ 
fore  to  the  lowest  order  we  can  approximate  af  ~  aj, 
and  (1)  reduces  to  an  integrable  Hamiltonian  Hq  made 
of  non-interacting  bosons 

Hq  =  ^  Jij(a\aj  +  a-aj-  +  h.c.)  +  2B  ^  dja,,  (4) 

i<j  i 

H,  =  H  -  H0.  (5) 

Here  Hi  contains  interactions  between  the  bosons  which 
are  parametrically  small  in  n,  and,  as  a  result,  we  can 
treat  Hi  as  a  perturbation  to  Hq.  Thus,  it  is  natural 
to  expect  the  system  (in  the  thermodynamic  limit)  to 
first  relax  to  a  prethermal  state  described  by  the  GGE 
of  H0,  and  to  later  relax  to  a  thermal  state  described 
by  the  full  H.  Naively,  we  expect  the  thermalization  to 
happen  at  a  time  scale  much  longer  than  the  relaxation 
to  GGE,  based  on  the  different  energy  scales  of  H0  and 
Hi.  However,  this  is  not  always  the  case,  as  discussed 
below. 

To  explicitly  define  the  GGE  of  H0 ,  we  would  need 
to  first  diagonalize  Hq  and  find  the  integrals  of  motion. 


Note  that  (|4|  only  involves  Jij  for  i  <  j.  For  convenience, 
we  will  define  a  new  matrix  J  such  that  Ja  =  0  and 
Jij  =  Jji  =  Jij  for  i  <  j.  Hq  can  be  rewritten  as 


Hq  —  J ij 


a|a,- 


l(atafj  +  aw)  +2 


a\az. 


(6) 

An  orthogonal  matrix  V  can  be  used  to  diagonalize  the 
matrix  J  as  J2ij  VikJijVjk >  =  Vkhw ,  where  {vk}  are  the 
eigenvalues  of  matrix  J .  Introducing  ck  =  JT  Vikdi,  we 
have 


Ho  =  £ 


{vk  +  2  B)c\ck  +  -vk(c\c\  +  CfeCfe) 


(7) 


Next,  we  perform  a  Bogoliubov  transformation  ck  = 
cosh (9k)dk  —  sinh(0fc)dj(  with  0k  =  |  tanh-1)^  ^2b) 
fully  diagonalize  H0 : 

Hi)  =  ^  ekd\dk ,  ek  =  2 y/B(B  +  vk).  (8) 

k 


Since  integrable  models  have  an  extensive  number  of 
conserved  quantities  that  are  not  taken  into  account  by 
canonical  ensembles  from  statistical  mechanics  the  GGE 
was  developed  to  make  predictions  about  equilibrium  val¬ 
ues  of  observables  in  these  systems  by  incorporating  the 
additional  integrals  of  motion  [121 H31  [33-SO] .  The  GGE 
for  Hq  is  defined  as 


Pgge  = 


g  12k  ^ kdj.dk 

Tr(e_2d'=  Afc4  dfc)’ 


(9) 


with  AfcS  determined  by  (d\.dk)o  =  ( d\dk)cGE ,  where 
the  notation  (•  •  •  )o  denotes  the  expectation  value  in  the 
initial  state  1^0)5  and  the  notation  (•  •  •  )gge  denotes  the 
expectation  value  in  pgge ■  We  observe  relaxation  to  the 
GGE  if  for  any  local  observable  O ,  ( 0{t ))  «  tr (Opgge)- 
The  GGE  is  equivalent  to  the  grand  canonical  ensemble 
for  a  general  quantum  system  when  the  only  conserved 
quantities  are  particle  number  and  total  energvjT2l. 

Using  the  fact  that  our  initial  state  \i/jq )  is  always  a 
Fock  state  in  the  basis  of  {a-(q},  the  value  of  (d\.dk)o 
can  be  calculated  by  the  following  formula 


(dldk) 0  =  cosh(26»fc)  ^  Vl2fe(ajal)0  +  sinh2(6»fc).  (10) 

i 


To  calculate  the  expectation  values  of  of  =  2 aja^  —  1 
in  the  GGE,  we  use  the  following  equations 

{a\  ai)cGE  =  <£vifcV^4cfe,>GG£ 

k,k' 

=  (cosh  {29k)dkdk  +  sinh  2{0k))GGE 

k 

=  54  Vik[cosil(26k){d{dk) 0  +  sinh2(0fc)], 

k 


(11) 


where  we  use  the  fact  that  pgge  is  diagonal  in  the  Fock 
basis  of  {44},  so  (d\dk>)GGE  =  (44')ggb  =  0  for 


k  ±  k'. 


where  we  have  set  w2  =  4,^7  as  an  “effective”  axial 
trapping  frequency.  The  dimensionless  matrix  K  charac¬ 
terizes  the  dipolar  interactions  between  ions: 


Single-particle  properties  of  Ho 

Since  vk/B  is  small,  we  can  expand  ek  and  4  in  vk/B 
to  the  leading  order: 

efe  «  2 B  +  uk,  dk  «  ^Vifc(al  +  £^aj).  (12) 

i 

This  means  to  understand  the  single-particle  properties 
of  H0  (Eq.|§.  we  just  need  to  understand  the  properties 
of  the  eigenvalues  {vk}  and  eigenvectors  V  of  the  matrix 
J .  We  emphasize  that  the  matrix  J  defined  in  (|3])  differs 
from  the  matrix  J  used  in  defining  H0,  because  Ju  7^  0 
by  the  above  definition  ([3]).  This  is  because  has  no 
physical  consequence  in  the  Ising  Hamiltonian  as  (erf)2  = 
1. 


Eigenvalues  and  eigenvectors  of  J 

Let  us  first  try  to  understand  the  properties  of  the 
eigenvalues  and  eigenvectors  of  the  matrix  J.  To  make 
this  possible,  we  will  need  to  approximate  the  spacing 
between  ions  to  be  uniform.  While  this  is  not  true  in  the 
current  experiment  due  to  the  harmonic  trapping  poten¬ 
tial,  the  inhomogeneity  in  ion  spacing  is  not  responsible 
for  the  observed  prethermalization  |14|  and  from  now  on 
we  will  assume  that  the  ions  are  equally  spaced. 

We  now  write  down  the  motional  Hamiltonian  of  N 
ions  trapped  along  the  2:  direction  ignoring  the  ions’  mo¬ 
tion  along  the  y  direction  for  simplicity  since  we  barely 
couple  the  spins  to  the  phonons  in  that  direction: 


Hm=  Ef= 


Pi,X 
2  M 


Pj,z 

2  M 


Q 

47re0 


+  V{zi)  +  \Mu2xx2 

EN  i —  1 

2=1  Z— /  7  =  1 


1  yj  (Zi—Zj  )2+(xi—Xj 


(13) 

(14) 


Here  {xi,  Zi,pitX,pitZ}  are  respectively  the  coordinates 
and  momenta  of  the  ith  ion  in  the  x  and  2  directions.  M 
and  Q  are  the  mass  and  charge  of  each  ion,  and  u jx  is  the 
transverse  trapping  frequency.  The  ions  will  be  equally 
spaced  with  a  spacing  ao  if  V(z)  =  —  47Joaii  log(l  — 
22/L2),  with  L  =  lVao/2. 

By  expanding  the  Coulomb  interaction  around  the 
ions’  equilibrium  positions  up  to  second  order  in  posi¬ 
tion,  the  motional  Hamiltonian  in  the  x  direction  can  be 
written  as: 


//,, 


Pi,x 

^  2 M 


N 


-M 


Lolxf 


ujz  KijXiXj 

i,j= 1 


(15) 


K*j  =  -\  i~j\-3,  Kii  =  -YJKij,  (16) 

jji 

The  exact  analytical  expressions  for  the  eigenvalues 
{Km}  and  eigenvectors  {V),m}  of  K  cannot  be  obtained. 
But  we  can  employ  a  first-order  perturbation  theory  and 
assume  that  the  eigenvectors  of  K  are  approximately  the 
same  as  those  of  a  nearest-neighbor  coupling  matrix  K. 
As  a  result: 


f  4W  171  = 

cos[^P(i  —  §)],  m  =  l,2,  •••TV— 1, 

^4  2  —  2  cos  wJtl  .  .  .  . 

- r3  ,  (m  =  0, 1,  ■  ■  ■  Af  —  1).  (18) 

r= 1 

Note  that  the  (*  —  i)  above  ensures  that  the  phonon 
modes  are  either  symmetric  or  antisymmetric  under  the 
spatial  inversion  of  the  chain  (i  — >  N  +  1  —  i). 

As  a  result,  the  eigenvectors  of  the  matrix  J  are  given 
by  and  the  eigenvalues  are  given  by 


A 


m 


h(Sk)2n2 

2 M{p2  -  w2 +w2kto)’ 


(19) 


Finally,  we  point  out  that  importantly,  Ju  is  in  general 
non-uniform.  This  can  be  seen  in  the  following  two  limits: 

1.  When  p2  —  w2  co2tim  for  all  m,  we  expect  to 
decay  as  1/| i  —  j |3  (a  — >  3  limit),  and 


h(Sk)2td2u}2 
2M (p2  -  Jxy 


(20) 


Ju  in  this  limit  is  very  close  to  uniform  in  the  large 
N  limit,  except  for  i  close  to  1  and  N  (see  Fig.  i>)- 

2.  When  p2  —  oj2  <C  w2Km  for  all  m  7^  0  (a  — »  0 
limit),  we  can  separate  out  the  m  =  0  term  and 
approximate  Ju  by 


2  Mix2  N\  +  ^ 


(21) 


Note  that  even  in  the  large  N  limit,  Ju  is  non- 
uniform  across  the  entire  ion  chain  (see  Fig.  i>)- 
An  analytical  formula  can  be  obtained  if  we  ap¬ 
proximate  Km  by  including  only  the  r  =  1 

(nearest-neighbor)  term  in  Eq.  18  leading  to  Ju  ~ 
~  ^)2  +  constant. 
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Eigenvalues  and  eigenvectors  of  J 


As  shown  above,  when  the  interactions  described  by 
Jij  are  very  long-ranged,  can  be  rather  non-uniform. 
This  will  result  in  a  qualitatively  different  structure  be¬ 
tween  the  eigenvalues  and  eigenvectors  of  the  two  ma¬ 
trices  J  and  J .  To  see  this,  we  first  notice  that  the 
eigenvalues  and  eigenvectors  of  J  are  similar  to  those  of 
the  Hamiltonian  for  a  free-particle  in  an  square  well  po¬ 
tential.  This  connection  can  be  formalized  by  going  into 
the  continuum  limit  and  introducing  a  continuous  mo¬ 
mentum  q  =  mir/N  £  (0, 7r ) .  The  eigenspectrum  A (q)  of 
J  is  minimized  at  q  =  7r.  Expanding  A (q)  around  q  =  ir 
and  using  Eq.[T8fT9}  we  obtain 


A (q)  ~  0({q  -  tt)°)  +  ^P~(q  ~  tt)2  +  0[(q  -  tt)4], 


Meff  =  M 


\li2  —  to2 
—  ji  r  Lf  & 


2  Meff 


■4C(3)]‘ 


to2n2  In  2 


(22) 

(23) 


corresponding  to  the  dispersion  relation  of  a  massive  par¬ 
ticle  with  an  effective  mass  Me f  f  and  an  effective  momen¬ 
tum  ( Sk)q . 

The  Schrodinger  equation  for  the  above  particle  can  be 
written  in  the  position  space  parameterized  by  a  contin¬ 
uous  coordinate  z  £  [1,  iV]  that  replaces  the  discrete  ion 
index  i  £  {1,  2,  •  •  •  Nj: 


h2(Sk)2  d 2 

2 Meg  dz2 


i’(z)  =  Ei/j(z), 


(24) 


with  the  boundary  condition  %l>(z  <  1)  =  tl>(z  >  N)  =  0 
corresponding  to  that  of  a  particle  in  a  box  potential. 
Here  the  eigenwavefunction  ipm(z)  ~  (— 1  )lVitTn  (for  z  = 
i),  and  the  eigenenergy  Em  «  h\m  (up  to  a  constant 
shift)  near  q  =  tt. 

We  can  similarly  map  the  eigenvalue  equation  of  J  to 
a  Schrodinger  equation  of  a  massive  particle: 


-w£'s(2)+t,w=^w-  (25) 

However,  we  now  have  an  effective  potential  U{z)  due  to 
the  fact  that  Ju  =  0  ^  Ju,  and  up  to  a  constant  energy 
shift 


U(z) 


oo  z  <  1  or  z  >  N 
-  Ju  z  =  i  £  {1, 2,  -  -  •  ,  N} 


(26) 


As  discussed  in  the  previous  section,  for  fj,2  —  uj2  to2Km 
(where  a  — >  3),  the  potential  U(z)  will  be  nearly  flat,  and 
the  eigenvalues  and  eigenvectors  of  J  are  similar  to  those 
of  a  particle  in  a  box.  However,  for  p t 2  —  to2  <C  to2Km 
(where  a  — ►  0),  U{z)  has  the  shape  of  a  double  well 
potential  (Extended  Data  Fig.  io- 


i/N 


Extended  Data  Figure  3.  (a)  The  diagonal  matrix  element 
—  Ju  [in  arbitrary  units  and  shifted  by  min(Jjj)]  that  deter¬ 
mines  the  single-particle  potential  U{z)  for  N  =  100  ions.  We 
choose  the  parameters  that  make  Jij  decay  approximately  as 
l/r“  with  different  values  of  a  shown  in  the  plot.  As  a  de¬ 
creases,  the  potential  changes  continuously  from  nearly  flat 
to  an  approximately  harmonic  anti-trap.  Together  with  two 
hard  wall  potentials  at  i/N  =  0, 1,  the  potential  looks  like  a 
double  well  that  becomes  deeper  for  smaller  a.  (b)  The  eigen¬ 
vector  Vj.io  corresponding  to  the  10th  lowest  eigenvalue  of  J 
for  a  N  =  100  ion  chain.  For  small  q’s,  the  eigenvector,  as 
well  as  the  wavefunction  'F(z)  of  Eq.  (25 1,  is  localized  inside 
the  two  wells.  For  large  a’s,  the  eigenvector  is  delocalized  and 
similar  to  that  of  a  particle  in  a  box. 


Discussion 


Ignoring  tunneling  between  the  two  deep  wells  of  a 
double-well-shaped  potential,  the  low-energy  eigenstates 
of  a  massive  particle  in  such  a  potential  are  localized 
orbitals  inside  either  well  (Extended  Data  Fig.  !>)■  The 
tunneling  rate  (which  is  exponentially  small  in  the  height 
of  the  barrier)  splits  the  degeneracy  of  the  localized  or¬ 
bitals  in  each  well  and  leads  to  pairs  of  symmetric  and 
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antisymmetric  (upon  the  spatial  inversion  of  the  chain) 
wavefunctions.  This  physical  picture  explains  the  ob¬ 
served  prethermalization:  Initial  excitations  placed  in  the 
left  half  of  the  chain  will  be  localized  for  an  extended 
period  of  time  under  the  evolution  of  Ho ,  until  the  tun¬ 
neling  between  the  two  wells  eventually  delocalizes  the 
excitations. 


Importantly,  the  double- well-shaped  potential  is  emer¬ 
gent  in  our  system,  because  the  non-interacting  Hamil¬ 
tonian  Ho  does  not  contain  an  inhomogeneous  poten¬ 
tial.  This  emergent  inhomogeneity  is  somewhat  a  sur¬ 
prising  effect  because  the  motional  Hamiltonian  (Eq.  14 ) 
and  spin-motion  couplings  induced  by  the  Raman  lasers 
are  all  homogeneous.  The  key  reason  is  that  the  long- 
range  interactions  break  the  translational  invariance  of 
the  Ising  Hamiltonian  (even  in  the  thermodynamic  limit). 
This  is  in  contrast  to  an  open-boundary  spin  chain  with 
short-range  interactions,  where  it  is  safe  to  assume  trans¬ 
lational  invariance  for  sufficiently  large  system  sizes. 


The  notion  of  a  boundary  starts  to  break  down  for  suf¬ 
ficiently  long-ranged  interactions,  and  therefore  we  can¬ 
not  attribute  the  observed  prethermalization  to  bound¬ 
ary  effects.  Usually,  boundary  effects  only  affect  a  finite 


number  of  eigenstates  and  do  not  affect  local  quenches 
in  the  bulk.  However,  here  there  is  extensive  number  of 
eigenstates  that  are  localized  in  the  two  wells  described 
by  the  potential  U(z),  and  excitations  placed  an  exten¬ 
sive  number  of  lattice  sites  away  from  the  edges  are  still 
subject  to  qualitatively  similar  dynamics. 

Finally,  we  point  out  that  interactions  in  Hi  can  also 
delocalize  the  initial  spin  excitations  placed  in  one  of  the 
wells,  and  eventually  thermalize  the  system.  As  a  result, 
there  is  an  interesting  interplay  between  the  timescales 
of  prethermalization  to  the  GGE  and  of  thermalization. 
If  the  interactions  in  Hi  are  sufficiently  weaker  than  the 
kinetic  tunneling  rate  in  Hq,  which  can  be  achieved  by 
increasing  the  magnetic  field  strength  B  or  changing  the 
range  of  interactions,  then  we  expect  the  system  to  have 
two  prethermal  phases  before  thermalization,  with  the 
observed  prethermalization  followed  by  prethermaliza¬ 
tion  to  the  GGE  of  Hq.  If  instead  the  tunneling  rate 
is  sufficiently  smaller  than  the  interactions  in  H i,  then 
the  prethermal  states  described  by  the  GGE  of  Ho  may 
never  appear  during  the  time  evolution.  These  interest¬ 
ing  multi-stage  relaxation  processes  will  require  future 
experimental  investigations  with  longer  coherence  times 
and  larger  spin  chains. 


